
alpha = -180:5:180;
beta = -90:5:90;

alpha_r = deg2rad(alpha);
beta_r = deg2rad(beta);

[alpha_mesh, beta_mesh] = meshgrid(alpha_r, beta_r);

%R = abs(cos(alpha_mesh).*1);
R = [1/4 -1/4 0;
     -1/4 -1/4 0;
     -1/4 1/4 0;
     1/4 1/4 0];
 alpha_f = 0;
 beta_f = deg2rad(30);
 K = [0, cos(beta_f), sin(beta_f)];
 dPhi = 2*pi*K*R';
 B = exp(-1*dPhi);
 for i=1:size(alpha_mesh,1)
     for j=1:size(alpha_mesh,2)
         alpha = alpha_mesh(i,j);
         beta = beta_mesh(i,j);
         K = [cos(alpha)*sin(beta), sin(alpha)*sin(beta), cos(beta)];
         B1 = exp(1*2*pi*K*R');
         DN(i,j) = abs(sum(B.*B1)^2);
     end
 end
 
 
[x, y, z] = sph2cart(beta_mesh, alpha_mesh, DN);

figure(1);
surf( x, y, z)
xlabel('x'); ylabel('y'); zlabel('z')

minc = min( [min(min(x)) min(min(y)) min(min(z))] );
maxc = max( [max(max(x)) max(max(y)) max(max(z))] );
xlim([minc maxc]); ylim([minc maxc]); zlim([minc maxc]);

% surf( alpha_mesh, beta_mesh, R)